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1. Introduction 

The interaction of strong gravitational and magnetic fields is important in a variety of 
astrophysical phenomena. The Blandford-Znajek process | Ijj, in which a magnetized plasma 
can extract spin energy from a black hole, is a promising mechanism for understanding 
relativistic jets in AGNs, galactic microquasars, and gamma-ray bursts. Furthermore, the 
magnetorotational instability is an important mechanism to transport angular momentum in 
accretion disks 0. Neutron stars and pulsars may also have intense magnetic fields, with 
magnetars being an extreme example. Magnetars are models for soft gamma-ray repeaters 
and anomalous x-ray pulsars. 

A large body of work in relativistic magneto-hydrodynamics has been done in flat space 
or special relativity. The incorporation of significant gravitational effects, however, requires 
general relativity. The Arnowitt-Deser-Misner (ADM) formulation of the Einstein equations 
is adapted to solving the initial value problem in general relativity 0, and pioneering work 
in GRMHD using the ADM formulation was done by Sloan and Smarr 0, and Evans and 
Hawley 0. A renewed interest in GRMHD is evident from recent work on both fixed 
backgrounds ID[Zl[8|0IIo|IID and with dynamic geometries 02] 021]- 

In this paper we describe our method for solving the relativistic MHD equations. The 
equations are derived for a general, arbitrary spacetime. The standard flat space test problems 
are performed on multiple grids with an excised region. We investigate hyperbolic and elliptic 
divergence cleaning for the MHD equations. Although the numerical results shown here are 
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done in flat space, we develop techniques applicable for black hole spacetimes with excision. 
Secondly, we have designed our method for seamless integration with the Einstein equations 
when using adaptive mesh refinement (AMR). The combination of general relativity and 
MHD in evolutions with AMR will be presented in subsequent papers. 

The fluid equations are fundamentally conservation equations, and this conservation 
property can be expressed numerically by using finite volume (FV) discretizations. Here the 
domain is discretized into volume elements or cells of finite size, and the evolution variables 
represent cell averages, such as energy or momentum. Essentially Non-Oscillatory (ENO) 
numerical schemes however, allow either a FV, or, for uniform grids, a finite difference 
(FD) discretization. FD functions represent point values at discrete points as opposed to the 
averages in FV schemes. A FD formulation of the MHD equations in general relativity is 
compelling for two reasons: (1) communications between arbitrary grids are simplified as 
only point values are required; and (2) solving the combined fluid and Einstein equations 
with AMR is simplified when both sets of equations are discretized in the same manner. The 
simultaneous refinement of both vertex centered (FD) and cell centered (FV) grids results in 
staggered grids. As the Einstein equations are frequently discretized with finite differences, 
a combined general relativity and MHD code can be simplified if both sets of equations are 
discretized in the same manner. We choose a central Convex ENO (CENO) method with FD 
discretization for the MHD equations. 

Black hole excision is commonly used in numerical studies of black hole spacetimes. 
The inner excision boundary must be carefully constructed such that numerical modes do not 
enter the domain through the inner boundary 115111611771 . Coordinate systems adapted to the 
horizon’s geometry allow one to excise the largest volume of spacetime. This is advantageous 
because gradients in the gravitational fields become larger the closer one is to the singularity, 
requiring greater computational resources to adequately resolve them. Finding a global 
coordinate system adapted to all boundaries may only be possible for the most symmetric 
cases. Thus, multiple coordinate patches may be necessary to cover the entire domain. 

Different coordinate systems can be implemented computationally using multiple grids 
with appropriate communication defined between grids. One approach is to use touching or 
abutting grids, in which all boundary points on neighboring grids coincide. These grids have 
been used in black hole evolutions with both spectral ca and finite difference methods ED. 
A second approach uses multiple grids that overlap j20l BTS . Information between grids is 
communicated via interpolation. Overlapping grids allow for greater freedom in choosing 
numerical schemes, coordinate systems, and moving some grids with respect to others. The 
feasibility of using overlapping grids for moving black holes was explored by successfully 
solving the Klein-Gordon equation on fixed black hole backgrounds for both highly boosted 
(v = 0.98c) Schwarzschild C3 and Kerr black (a = 0.99) holes t22l . It is most natural to 
implement HRSC schemes for fluids on overlapping grids. 

The magnetic field B is evolved with the MHD fluid equations and must also satisfy 
the constraint V • B = 0. Experience has shown that error in this constraint will grow to 
unacceptable levels (see, for example, figure 0 and discussion below), leading to unphysical 
solutions, unless the constraint is actively enforced. Moreover, the constraint growth is 
exacerbated when weakly hyperbolic formulations of the MHD equations are used ED 
Some techniques to enforce this constraint include constrained transport, projection methods, 
and hyperbolic divergence cleaning t24ll25)f Constrained transport differences the evolution 
equations for B such that V • B = 0 is satisfied to machine precision for one particular 
discrete divergence operator. Naturally, the continuum constraint is only satisfied to the level 
of truncation error, which can be easily seen by evaluating the divergence with a different, 
consistent discrete divergence operator |24j. While some constrained transport schemes may 
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be used with structured AMR 1261 , we turn to other methods for greater freedom in choosing 
multiple coordinate systems. Hyperbolic divergence cleaning adds a new field designed such 
that divergence errors are propagated off the grid ED, and is similar to the A-system of 
Brodbeck et al for the Einstein equations |28j. Projection methods involve solving an elliptic 
equation for a correction to B, such that it satisfies the constraint. Both the projection method 
and hyperbolic divergence cleaning are easily implemented with overlapping grids. 

This paper presents details of our method and gives results of numerical tests. The MHD 
equations in general relativity are derived in section |2] Section 0 presents the numerical 
scheme. Section 0 discusses divergence cleaning for the MHD equations. All numerical 
results are performed on overlapping grids, and are presented in Section^] 

2. The MHD equations in general relativity 

We first derive the equations of motion for relativistic MHD and a dynamic spacetime. The 
equations are written in conservation form as required for High-Resolution Shock-Capturing 
(HRSC) numerical methods. We then discuss the transformation between conserved and 
primitive variables. 

2.1. Equations of motion 

To begin, we assume a stress energy tensor of the form 

Tab = [po (1 + e) + P] U a llb + Pg a b + FacPb — ^.9 a b F c dF Cd , (1) 

where the first few terms describe the fluid and the final two terms the electromagnetic field. 
The fluid and electromagnetic components are coupled through the relativistic form of Ohm’s 
law: 

J a + (u b J b ) u a = aF ab u b , (2) 

where J a is the 4-current. The ideal MHD approximation is simply the statement that the 
fluid has perfect conductivity, i.e., a —> oo. Equivalently, this can be expressed as 

F ab u b = 0, (3) 

which states that the electric field in the frame of the fluid vanishes. This is sometimes referred 
to as the “freezing-in” condition of the magnetic field; namely, in the frame of the fluid, the 
magnetic field lines are frozen to the fluid and carried along with it. 

With this in mind, a convenient set of substitutions for the electromagnetic variables is 
to define 4-covariant “electric” and “magnetic” four-vectors 

e a = F ab u b , b a = *F ab u b , (4) 

where *F ab = e abcd F c d /2 and e abcd is the standard totally antisymmetric Levi-Civita tensor. 
Note that we can write these as 

Fed. tl c Cd C-cdeftt b ^, * Fed — tlcbd Ccdeftl 6^,(5) 

where we have the constraints u a e a = 0 = u a b a . All the information in the Maxwell tensor, 
F a b , is now contained in these two four vectors. 

With these substitutions, the electromagnetic part of the stress tensor can be written as 

T ab ti a Ub [e c e -j- b c b ] -f- — g ab [ec^ T b c b ] e-a^b b a b b T 2 'u^ a 65 ^ c ^ e e u b . (6) 
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In the MHD approximation, the electric four vector is identically zero and the full stress tensor 
for MHD can be written as 

Tab = [po (1 + e) + P + b c b c ] Ua'Ub + P + — b c b c g a b — b a bb . (7) 

The matter equations of motion can now be written in conservation form 

V a T ab = 0, V a *F ab = 0. (8) 

To these must be appended the baryon conservation equation V a ( pou a ) = 0. 

In a general spacetime we decompose these equations in the usual ADM 3+1 split by 
projecting along and orthogonal to a unit normal vector, n a , which is orthogonal to a foliation 
of spatial hypersurfaces. The projection tensor is 

hab = Pab T tlaTlb , (9) 

with g a 6 the metric on the 4-manifold. The Einstein equations have the usual 3+1 form with 
both evolution and constraint equations. Because our focus in this paper is developing a 
robust MHD code, we will emphasize and solve the flat spacetime equations in later sections. 


However, our approach in deriving the equations in this section is completely general. 
Conservative variables are defined in the conventional way 

E =T ab n a n b , (10) 

S b = - T ac n a h b c , (11) 

(J -T) cd =T ab h a c h b d . (12) 

With respect to the MHD stress tensor, these give 

E = [p 0 (l + e) + P + b c b c ](n a u a ) 2 - P +\b c b c - {n a b a f , (13) 

S b = - \p 0 (1 + e) + P + b c b c ] (n a u a ) (±u) b + (n a b a ) {±b) b , (14) 

(±T) cd = [p 0 (l + e)+P+b c b c ](±u) c (±u) d + P+\b c b c h cd -(±b) c (±b) d , (15) 
where we have defined 

W = -n a u a , v a = - (_L u) a , (16) 


and (i. X )" = h a bX b denotes a projection. Note that W is the Lorentz factor between the 
fluid frame and the fiducial observers moving orthogonally to the spatial hypersurfaces. In 
addition, v a is the (purely spatial) coordinate velocity of the fluid. The matter equations are 
projected along and orthogonal to n a , and expressed in terms of the conserved variables 

0 = - n a 8 a E + KE - \D a (a 2 S a ) + (lT) at K ab , (17) 

a z 

0 = hbc \-n a d a S b + I\S b + 2S a K a b - -S a d a (3 b - -D a (a (±T) ab ) - — e] , (18) 

0 = D a ( *F ab n b ) , (19) 

0 = h bc ~n a d a (*F de n d h b e ) + *F db n d I< + -D a (a (J_ *F) ab ) - - *F da n d d a f3 b , (20) 

a a 

0 = -n a d a ( aD ) + - D a ( aDv a ) - KW, (21) 

a a 

where a and 3 b are the ADM (3+1) lapse and shift, K ab is the extrinsic curvature, and D a 
is the covariant derivative compatible with h ab - These equations, in order, are the energy 
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equation, the Euler equation, the no monopole constraint, the induction (or Faraday) equation 
and the baryon conservation equation. 

It is advantageous to use the standard magnetic field as the evolution variable, rather than 
the magnetic four vector b a . This amounts to working in the frame of the fiducial observers 
moving along n a instead of in the fluid frame. The electric and magnetic fields in this frame 
are then 

E a = h a »F bc n<, B a = \e abc F^. (22) 

where e a bc = n d Cdabc■ The ideal MHD approximation then becomes a relation giving the 
electric field in terms of the magnetic field in the frame of the orthogonally moving observers: 

Ea = — e abc u b B c . (23) 

n d u a 

In practice, two modifications are made to the MHD equations in order to solve them. 
First, we evolve the quantity r = E — D instead of E alone. This is often done to have an 
energy quantity that reduces to the Newtonian value in the nonrelativistic limit. Secondly, the 
source term in the induction equation can be eliminated by combining that equation with the 
no-monopole constraint. The final form for our matter equations thus becomes 


d t yv h rj + di 

d t (VhSb) + d, 


V~9 [ S i ~ — t — v i D 




= V=9 


(±T) ab K a b — — S a d a a 
a 


sFV) (±T)' b -Z-S b 


= V=~9 


3 r; h (irr,+ -S a d b (3 a - —d b a E 


a 


a 


‘ a(^ B .)=o. 


Vh 

d t (VhB b ) +8i 




V=g B b i V * - z- - B i v b - 


= 0 , 


d t ( 'VhD ) + di V~gD ^ ^ 


= 0. 


(24) 

(25) 

(26) 

(27) 

(28) 


2.2. Primitive and consen’ed variables 

The evolution equations give the time dependence of the conserved variables, u = 
(D, Si,T, Bj) T , but they also depend on the primitive variables w = (po,Vi, P,bj) T ■ As 
discussed in this section, for relativistic fluids the transformation from conserved to primitive 
variables is transcendental. The ability to solve for physical values of the primitive variables 
under a wide variety of conditions is an important and challenging part of writing a relativistic 
fluid code. 

The conserved variables are 

D = Wp 0 , (29) 

S b = (h + b c b c ) W 2 v b + ( n a b a ) (_L b) b , (30) 

r = {h + b c b c ) W 2 -P~ hj c b c - (n a b a f - Wp 0 , (31) 

“ = -Wb a - u a ■ ( n c b c ), (32) 


B‘ 
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where the fluid enthalpy is h = po(l + e ) + P- To obtain the inverse transformation, we 
reduce the problem to the solution for the roots of a single nonlinear function. The method is 
as follows. 

We eliminate the magnetic four vector, 6\ from the above equation using 

6“ = --^ [B a +u a -(±u) b B b \ . (33) 

On replacing this, we get 


D = Wp 0 , 

Si= (hW 2 + B 2 ) Vi - (B j Vj ) Bi, 


r 


= hW 2 + B 2 - P- 


(- B l Vi y 


w 2 


Wpo, 


(34) 

(35) 

(36) 


where B 2 = BiB \ v 2 = VjV 1 , and the indices are raised and lowered by the spatial metric 
hjj. The spatial norm of v l can be expressed in terms of the Lorentz factor 


W 2 


1 

1 — V l Vi 


(37) 


Density and pressure, two primitive variables, can be expressed as 

Po = D = D s/l - v 2 , P=(h-p 0 )^Y~- (38) 

Note that we assume in this section a T-law equation of state. 

It now remains to find v l (or W) and h from our knowledge of D, Si, t and Bi. We 
contract B l with Si 


SiB i = hW 2 (B i Vi), 


(39) 


and use this to eliminate B l Vi in the expressions above for r and Si. From S' 1 Si we derive the 
expression 

- (hW 2 ) 2 W 2 SiS i + ( hW 2 ) 2 ( hW 2 + B 2 ) 2 (W 2 - l) 

- W 2 {2hW 2 + B 2 ) ( S i B i ) 2 = 0. (40) 

This can be solved for W 2 in terms of conservative variables and the quantity x = hW 2 : 


W 2 = 


1 - 


(2a; + B^jBiSj) 2 + x 2 (S->Sj) 


n -1 


(41) 


x 2 (x + B 2 ) 2 

Finally, we substitute <ED into the equation for r (which comes about on using our above 
expressions for the density and pressure): 

x 2 = i (BJSj) 2 .(42) 


x ( 1 — 


r-1 i 


r w 2 


-D 1 - 


r-i i 


r w 


i 


777 - r+ o B2 ( 1+v2 ) 


The full expression is thus a nonlinear function in x, the roots of which we must calculate. 
Note that all the coefficients in this expression are conservative variables that on numerical 
integration of the evolution equations will be known at a given time level. Once x is obtained 
by solving ED, it is then straightforward to find W 2 , v 2 , h, po. P and b a . Equation ED is 
solved for x numerically using a combined Newton-Raphson and bisection solver. In practice, 
a floor is placed on po and P, and a typical value for the floor is 10~ 10 . The code simply halts 
when the primitive variable solver fails, and we do not interpolate values of the primitive 
variables from neighboring points. 
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3. Numerical methods 

This section describes the numerical methods used to integrate the MHD equations. The 
fluid equations are solved using the Convex ENO (CENO) method |29 , 30j. This method 
is based on point values (FD discretization) rather than cell averages (FV discretization), 
simplifying communications between grids. FD fluid methods are advantageous for multiple 
domain problems in general relativity when the Einstein equations are discretized with finite 
differences. A general relativistic MHD code with AMR, for example, can be simplified if the 
fluid and geometric variables are refined in the same manner. 

A second advantage of the CENO scheme is that the extension to systems of equations 
uses a component-wise decomposition, rather than one based on characteristic fields. This 
eliminates the need to calculate left and right eigenvectors of a Jacobian matrix. Centered 
schemes are more diffusive than those based on characteristic decompositions, but are easier 
to implement numerically and more efficient. Recent results show that these methods work 
well with relativistic fluids 13II1321 . Furthermore, the spectral decomposition of the Jacobian 
matrix for MHD is complicated by the existence of various degeneracies, and centered 
schemes have been widely used in relativistic MHD |8ll33l[T2lll3l . 

In a free evolution of the MHD equations, the divergence of the magnetic field can 
become very large. Constrained transport, a discretization technique for the magnetic field 
equations, is sometimes used to enforce the V • B = 0 constraint for the MHD equations. 
As we must interpolate data between arbitrary grids, we investigate two alternative methods 
for controlling this error. The first method uses an additional hyperbolic field for divergence 
cleaning, and the second is an elliptic projection method. Some comparisons are made using 
both techniques for relativistic fluids. 

Finally, we solve the equations on overlapping grids to facilitate the use of uniform 
grids in complex geometries. While the tests presented here are done in flat space, we use 
overlapping grids that mimic those used for excising black holes. 

3.1. Overlapping Grids 

In many problems it is necessary or advantageous to choose coordinate systems adapted to 
the boundaries. Except for highly symmetric systems, it is often difficult to choose a single 
coordinate chart appropriate for the entire problem. 

A typical case of interest is a spacetime containing a black hole. When black 
hole excision is used to remove singularities from the computational domain, adapting 
the boundary to the geometry of the event horizon yields the maximal excision volume. 
Coordinates adapted to the event horizon or one black hole, however, are usually not 
appropriate for large domains or for binary black holes. Thus multiple domain numerical 
techniques are appropriate in numerical relativity using both touching grids 1181 (T9! and 
overlapping grids l; 34ll35lfl7ll22l(36 1. We concentrate here on overlapping grids. 

To test computational techniques for black hole spacetimes, we perform all of the 
standard flat space tests in this paper on overlapping grids with excision. Figure Q shows 
the basic grid configuration used in these tests. The grid Q\ is a base grid from which a square 
region is excised about the origin, imitating the excision used around a black hole. A second 
grid, f/ 2 - is placed over the excision region to give the usual simply connected domain required 
for the flat space tests. Q 2 is rotated by an arbitrary angle 6 with respect to Q\. We choose 
coordinates (x,y) for Q\ and coordinates (£, 17 ) for C/ 2 . and the coordinates are related by a 
simple rotation 


£ = 2 ; cost? + y sin 9, 


rj = —x sin 9 + y cos 9. 


(43) 
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Figure 1 . This figure illustrates the overlapping grid structure used in all of the numerical tests 
presented in this paper. A region about the origin is excised from the base grid, Q i. A second 
grid, t/ 2 , covers the excision region and is rotated with respect to Q \. Data for the excision 
boundary of Q\ are interpolated from S2, and outer boundary data for S2 are interpolated from 
Si. Conventional out-flow boundary conditions are applied to the outer boundaries of Si. 


During the evolution boundary data at the inner boundary of Q\ are obtained by 
interpolating the solution from Qi. Outer boundary data for Qi are similarly obtained by 
interpolating the solution on Q\. For clarity figure ^indicates these interpolated points in a 
single row at each boundary. However, our third order evolution scheme requires a seven- 
point stencil, and we actually interpolate a band of three points at each inter-grid boundary. 
Finally, interpolation zones between grids are not allowed to overlap. 

Conservative variables are interpolated at grid interfaces, although a nonconservative 
interpolation scheme is used. These variables are typically smoother than the corresponding 
primitive variables, and have smaller relative jumps near discontinuities. The primitive 
variables are then recalculated from the interpolated variables. Simple Lagrangian 
interpolation can lead to oscillatory results near discontinuities, frequently resulting in 
unphysical states in relativistic fluid dynamics. (This effect may be less pronounced 
when using structured grids, such for Berger-Oliger AMR.) Therefore we interpolate with 
WENO interpolation ED. which is designed for use with discontinuous functions. WENO 
interpolation is summarized in |Appendix B| 

As mentioned briefly above, our FD fluid scheme on overlapping grids with WENO 
interpolation is not conservative. Conservative systems are often thought to be necessary for 
obtaining the correct weak solutions to the fluid equations. The effect of non-conservative 
boundary interpolation on such systems has been examined by Tang and Zhou f38l and 
Sebastian and Shu ED. Tang and Zhou found that convergent weak solutions can be obtained 
using nonconservative interpolation at grid interfaces. Sebastian and Shu found conservation 
errors of second order at grid interfaces for smooth solutions and first order errors for solutions 
with discontinuities at the interface. The numerical results presented in this paper are a direct 
demonstration that we are able to obtain the correct weak solutions when overlapping grids 
are used. 
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3.2. CENO 


High-Resolution Shock-Capturing (HRSC) methods are designed to solve hyperbolic 
conservation laws of the form 

d t u + d k f k (u) = s(u), u(0,x l ) = u 0 (x l ), (44) 

where u is a state vector, / k are flux functions, and s contains source terms. The semi-discrete 
discretization in one dimension is 


dui 

dt 


fi+l/2 ~ fi-1/2 
Ax 


s(u 


(45) 


where / is a consistent numerical flux. The accuracy of the scheme depends, 
among other things, on the estimation, or reconstruction, of u. i+[ / 2 from the quantities 
Ui-k, ■ ■■ .u,, ■ ■■. Ui +m , where k and m are integers. TVD schemes are named for the Total 
Variation Diminishing condition on the interpolating polynomial. These schemes accurately 
capture the dynamics of strong shocks without oscillations or Gibbs overshoot effects at 
discontinuities. The TVD condition, however, can also be overly restrictive, reducing the 
order of accuracy even at smooth extrema. 

The Essentially Non-Oscillatory (ENO) philosophy is that the interpolation stencil is 
chosen based on the local smoothness of the function. All points are used for smooth 
functions, and thus accuracy is not lost at smooth extrema. Near discontinuities the stencil 
is locally adjusted to use points away from the discontinuity. In relaxing the TVD condition, 
small oscillations may develop near discontinuities, but they are of 0(Ax k ), where k is the 
order of accuracy. The flexibility of the ENO philosophy has resulted in many extensions, 
including both FV and FD formulations, the Weighted ENO (WENO) approach, and schemes 
formally of very high order, etc. Shu reviews the different ENO methods in G3- 

The CENO method of Liu and Osher | [29j uses a FD discretization and Lax-Friedrichs 
flux splitting, eliminating the need for a characteristic decomposition. The interpolating 
polynomial in this scheme is produced by a convex combination of lower order interpolations. 
Near discontinuities this scheme is designed to produce results similar to TVD methods. This 
method was modified by Del Zanna and Bucciantini for relativistic fluids 1301 , and we follow 
their approach. 

We considered three different central or central-upwind numerical fluxes: (1) the Lax- 
Friedrichs (LF) flux, (2) the local Lax-Friedrichs flux (LLF) and (3) the Harten-Lax-van Leer 
(HLL) flux. We are primarily interested in highly relativistic systems, where the characteristic 
speeds approach the speed of light. In this limit the HLL and LLF fluxes reduce to the simple 
LF flux. Indeed, we have found very little difference between solutions calculated with the 
LF flux and those calculated with the LLF and HLL fluxes. Similarly, Del Zanna et al also 
reported nearly identical results from the LLF and HLL fluxes in this regime (33j. In this 
paper we use only the LF flux 


/ 


LF _ 
i+1/2 — 


f( U i+ 1/2) + f( U i+ 1/2) ( U i+ 1/2 U i+ 1/2) J 


(46) 


where itf +1/ / 2 and are the reconstructed states to the left and right of the interface 

at Si+ 1 / 2 , respectively. The CENO reconstruction for wfj _ 1/2 and /2 is described in 
|Appendix A| 

In FD ENO fluxes calculated from the point-wise values, f i+ i/ 2 , must be converted into 
consistent numerical fluxes, f i+1 / 2 . To order 0 (Ax fe ), the conversion from point valued 
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fluxes to conservative fluxes is given by 03® 

(k- 1)/2 


fi+ 1/2 = fi+x /2 + 51 a 2j (Ax) 2j (y) 


(47) 


i+l/2 


where a 2 = —1/24 and a 4 = 7/5760. For second order schemes the two fluxes are identical, 
/i+1/2 = fi + 1/2- A third order scheme, however, requires the correction 

fi+l/2 = (1 - Y^') fi+ 1/2’ ( 48 > 

where is a second-order non-oscillatory difference operator. I? 1 ' 2 //', is calculated, again, 
as a convex combination of the differences 

= minmod a°D { 0 2} fi, (49) 


where the minmod limiter is 

minmod(ai, a 2 , ■ ■ •) = 


minjafc} if all ak > 0 , 
maxjat} if all ak < 0 , 
0 otherwise. 


(50) 


and the one-sided (first-order) and centered (second-order) second derivative operators are 

D^fi = / i+2 -2/ i+1 +/ i , D^fi = fi — 2fi_i+fi_ 2 , D {2) fi = / i+1 -2/ i +/ i _ 1 .(51) 

The constants of are weights that may be chosen to bias towards centered differencing. Here 
we use or 1 = a 1 = 1 and a 0 = 0.7. 


3.2.1. Boundary conditions The third order CENO scheme has a seven point stencil. At 
physical boundaries three ghost zones are used, which are populated by simple extrapolation. 
For example, at a boundary x = xq we set it/ ■ = u 3 ■ for k = 0,1, 2 and all j. At inter¬ 
grid boundaries these ghost zones are set by interpolating from a covering grid using WENO 
interpolation, as described in |Appendix B| 

3.2.2. Time integration The semi-discrete equations (ED are integrated in time using third 
order Runge-Kutta. Various versions of RK3 exist, and we use the optimal third-order scheme 
of Shu and Osher that preserves the TVD condition 1391 

tx (1) =u n + AtL(u n ), 

u (2) = -u n + -u W + -AtL(u^), (52) 

4 4 4 

u n+1 = -u n + -u^ + -A tL(uW). 

3 3 3 v 2 

4. The V B = 0 constraint 

The magnetic field, B, must satisfy both its evolution equation as well as the solenoidal 
constraint: V • B = 0. Small numerical errors lead to violations of this constraint that, 
experience has shown, can grow rapidly. Left unchecked, violations of this constraint produce 
unphysical behavior. Various approaches can be used to enforce the solenoidal constraint, and 
here we consider two that can be applied to domains with multiple arbitrary grids and with 
a view to incorporating AMR: hyperbolic divergence cleaning and an elliptic projection of 
B. We do not consider constrained transport here because it requires that neighboring grids 
align in a structured manner, precluding its application to overlapping grids with arbitrary 
coordinates, resolutions and/or orientations. 
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4.1. Hyperbolic divergence cleaning 

Hyperbolic divergence cleaning is simple to implement numerically and comes with very little 
computational cost. A new scalar function ip is added to the system, which can be interpreted 
as a generalized Lagrange multiplier (GLM), and coupled to the magnetic field equations. 
The method can be implemented in various ways, and we follow the GLM method of Dedner 
et al for the classical MHD equations ED- See van Putten for another approach ED 

We specialize to Cartesian coordinates in flat space, and assume that %p satisfies a linear 
differential equation and couples to the magnetic field evolution equation according to 

d t B b + di ( B b v i - B i v b ) + gh-'djtp = 0, (53) 

Dtp + V • B =0, (54) 

where I? is a linear differential operator. Various choices for V can be made, giving 
hyperbolic, parabolic, and elliptic methods for divergence cleaning. We choose 

B>(tp) = \d t tp + \ ip, (55) 

< c l 

which combines some elements of both hyperbolic and parabolic operators, i.e., both 
propagating and damping the error. Equation then becomes 

dti!> + (? h V ■-& = -%!>. (56) 

c p 

Differentiating and combining shows that ip also satisfies the telegraph equation 

du-ip + C \d t ip - c\W 2 tP = 0. (57) 

C* 

Violations of the solenoidal constraint propagate with the speed Ch and the coefficient c p 
affects the damping rate. 

We have tested hyperbolic divergence cleaning using the cylindrical shock and relativistic 
rotor problems described below. In our tests we found that hyperbolic divergence cleaning to 
be very effective at keeping 11V • B| (2 bounded during an evolution. Some numerical tests are 
presented in section l5~4l We turn now to the elliptic projection method. 

4.2. The projection method 

Blackball and Barnes <42l first proposed an elliptic projection correction to the magnetic field 
such that it satisfies the constraint V • B = 0. The evolution equations are used to obtain a 
preliminary estimate for the magnetic field at the advanced time, B*. The corrected magnetic 
field at the advanced time, B n+1 , is then obtained by solving the system 

V 2 t/> = V B*, B ,i+1 = B* V'0. (58) 

Toth has shown for classical MHD that (1) V© is the minimal correction to B* that can 
be made such that V • B ra+1 = 0, and (2) the projection method gives the correct weak 
solution ED- 

The projection method can be implemented in different ways, e.g., the constraint can 
be imposed in Fourier space or physical space m We consider here only two different 
discretizations in physical space. Comparisons of different projection implementations, as 
well as to other techniques, such as constrained transport, are given by Toth |24|. and Balsara 
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and Kim E3. The divergence is discretized using the centered discrete operator D 0 , and the 
Laplacian can be discretized with either D 0 D 0 or D + D_. In one dimension the operators are 



Vi -(-1 Vi — i 

2 Ax 


(D 0 v) 


(59) 


The DqDq operator has a five point stencil in each direction, and the corrected magnetic field 
exactly satisfies the discrete divergence condition calculated using Dq. The D + D- operator 
has a three-point stencil in each direction, and the corrected magnetic field does not exactly 
satisfy any discrete divergence operator. As noted above, whether or not the magnetic field 
exactly satisfies a discrete divergence condition is something of a red herring: the continuum 
solution in both cases is only known to the level of truncation error. 

The projection method requires the solution of an elliptic equation dSU, which we solve 
iteratively after each complete Runge-Kutta cycle using a conjugate gradient or stabilized bi¬ 
conjugate gradient method. The error tolerance for the elliptic solver is set about 10 times 
less than the smallest expected truncation error, O (h 3 ). The solution of the elliptic equations 
is relatively efficient, requiring about 20-30% of the total run time for typical resolutions. 
Homogeneous Dirichlet boundary conditions are given for ip on all physical boundaries. At 
overlapping grid boundaries ip is interpolated as with other variables. 

When large corrections to B* are made in the projection process, problems can arise in 
reconstructing the corrected primitive variables. This may occur because both S a and r are 
themselves functions of the magnetic field, c. f. d36}. We compensate by also “correcting” 
these variables after projecting B, slightly modifying the conventional correction used in 
classical MHD 1241 . The projection algorithm can thus be summarized by 

(i) Solve the evolution equations for preliminary values at the advanced time u*; 

(ii) From u * calculate the primitive variables at the advanced time w/" + l («,*); 

(iii) Solve (ED for ip and compute B” +1 ; 

(iv) Recalculate u n+1 from w n+1 and B n+1 . 

This method assumes that the primitive variables more accurately reflect the correct solution 
in the projection process. Again, this is because p, v" , and P are not functions of the magnetic 
field, whereas the conserved variables S a and r are functions of B. 

5. Test Problems 

This section presents numerical results of our CENO scheme on overlapping grids. The 
first tests are a set of standard Riemann problem tests for relativistic MHD proposed by 
Komissarov 14311441 . Although these tests are inherently one dimensional problems, we run 
them on unaligned overlapping grids, making them effectively two dimensional problems. 
This allows us to test our divergence cleaning and interpolation methods on known solutions. 
We then discuss two test problems that are naturally two dimensional, the cylindrical shock 
and relativistic rotor problems described by Del Zanna et al (33j. Finally we present 
comparisons of the different divergence cleaning methods. 

5.1. Riemann problem tests 

Solutions of the Riemann problem are used to test shock capturing numerical methods, 
allowing one to verify that a method faithfully produces the fundamental shock, rarefaction 
and, for MHD, Alfven waves. The analytic solution of the relativistic MHD Riemann problem 
is given by Giacomazzo and Rezzolla EED. Komissarov presented several Riemann problem 
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Figure 2. This figure shows four Komissarov test problems. Solutions calculated on two- 
dimensional overlapping grids with elliptic divergence cleaning (plotted along the line y = 0) 
are compared with the one-dimensional solutions. The solutions are nearly indistinguishable. 
From left to right, top to bottom, these problems are: (1) Slow Shock, (2) Switch-off Fast 
Rarefaction, (3) Switch-on Slow Rarefaction, and (4) Compound Wave. Triangles indicate 
the solution on the excised base grid, Qi, and crosses indicate the solution on a grid covering 
the excision region, Q 2 , rotated 45° with respect to the base grid. Solid lines indicate the 
single-grid solution, and not the exact solution. 
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Figure 3. This figure shows two Komissarov test problems calculated on overlapping grids. 
The left frame shows Shock Tube # 1 and the right frame is Shock Tube # 2. See the caption 
to figure|3for further explanation. 


tests for RMHD, which we have used to test our code 14311441 . For comparison, results for 
Komissarov’s tests have also been published by other researchers 133111211X31 . With the initial 
discontinuity aligned with a coordinate direction, we are able to successfully reproduce all of 
Komissarov’s test problems, though we have reduced the Courant number to A = 0.4 in all 
tests, possibly because we are using the more diffusive Lax-Friedrichs flux. 

When the discontinuity is rotated with respect to the coordinates, we are able to run all 
Komissarov tests but two, the Fast Shock and the Collision. This problem occurs even when 
using a single computational grid, and thus is not related to using overlapping grids. In both 
cases we calculate unphysical values for the primitive variables, and the code is immediately 
halted before the solution is completed. We do not consider the Fast Shock problem here, 
and have modified the Collision problem by reducing the initial velocities from ±0.981 to 
±0.951. 

To provide the most comprehensive test of our algorithm, the Riemann problem tests 
presented here are performed on unaligned overlapping grids with excision. Excision is used 
here somewhat unconventionally, as the entire computational domain is simply connected. 
However, these Riemann problems can be used to test the excision and divergence cleaning 
algorithms. For example, we test that unphysical effects do not arise as waves pass through 
grid interfaces, and that the divergence cleaning methods do not adversely affect the solution. 
These solutions are then compared with those obtained from a single grid aligned to the initial 
discontinuity, with no excision and where divergence cleaning is unnecessary. Results from 
these tests are presented in Figures[2]and0 

In all tests F = 4/3, the Courant factor is A = 0.4, and elliptic divergence cleaning 
is used. Second order reconstruction with the minmod limiter is used in all tests except 
the collision problem, where first order reconstruction is more appropriate. The initial 
discontinuity in the fluid data is aligned with the coordinates (x, y) of a base grid, Q\. The 
region (x, y) £ [—4,4] is excised from Q\, and a second grid, Q 2 , covers this excision region 
to form the complete computational domain. Q 2 has coordinates (£, 77 ) £ [— i, -], which 
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are related to (x,y) by a rotation of angle 9 about the origin (I43> . We choose 9 = 45° for 
simplicity in plotting the results of our runs. In Figures [2] and [3] we plot data from the line 
y = 0 for Qi and the diagonal elements of Q 2 - 

In all cases we see that the solutions calculated on overlapping grids are nearly 
indistinguishable from those calculated on the single grid. In particular, we do not see 
reflections from the grid boundaries. Ill effects from elliptic divergence cleaning are also 
not observed in comparison to the single grid runs, where divergence cleaning is not used. 

5.2. Cylindrical blast wave 

Del Zanna et al presented this cylindrical shock test problem for relativistic MHD 1331 . The 
initial data consist of a uniform fluid background with p = 1, P = 0.01, v = (0, 0,0) and 
B = (4,0,0). Inside a disk of radius 0.16 centered at the origin, we set P = 1000. The 
adiabatic index is T = 4/3. An analytic solution is not available, but comparisons can be 
made with other published results f33lH~3l 

The solution is calculated using two overlapping grids with an excised region in the base 
grid. The base grid is uniform with (x, y) £ [—1,1], excluding the region (x, y) £ [—0.2,0.2]. 
The resolution is h = 0.008. A second grid covers the excision region with coordinates 
(£, rf) £ [—0.34, 0.34], rotated 50° with respect to the base grid. 

Figure |4] shows P and B 2 at t = 0.4. The solution is calculated using second order 
reconstruction, with a Courant number of A = 0.2, and elliptic divergence cleaning. The 
Courant number is lower than that used by Del Zanna et al , which may be a consequence 
of a different numerical flux and using a grid, Q 2 , rotated with respect to the initial magnetic 
field. The pressure difference is initially very large, leading to a strong out-going shock. 
The initially circular shock becomes distorted through interaction with the magnetic field, 
giving the elliptical profile observed in the figure. Comparing these results with those of other 
researchers, no artificial grid effects are observed in the solution, which could arise from using 
excision and overlapping grids. 

5.3. Relativistic rotor 

A second two-dimensional MHD test is the relativistic rotor, which evolves an initially rigidly 
rotating fluid in the presence of a magnetic field E3 Initial data consist of a constant 
background state with p = 1, P = 1, v = (0,0,0) and B = (1,0,0). Inside a disk of 
radius 0.1 at the center of the domain, the density is p = 10, and the fluid is rigidly rotated 
with ijj = 9.95. The linear velocity at the edge of the disk is 0.995 and W = 10. Finally, the 
adiabatic index is T = 5/3. 

The solution is calculated using two overlapping grids. The first grid has coordinates 
(x, y) £ [—1/2,1/2], excluding the region (, x , y) £ [—0.13, 0.13]. A second grid covers the 
excision region, and is rotated 27° with respect to the first grid. This grid has coordinates 
(£,77) £ [—0.2, 0.2]. The resolution on both grids is h = 0.0025. Figure |4] shows P and 
B 2 at t = 0.4. The Courant factor is A = 0.2, and second order reconstruction and elliptic 
divergence cleaning are used. Again, these results appear very similar to other published 
solutions, and no artificial grid effects are observed. 

5.4. Divergence Cleaning 

This section presents numerical results comparing the divergence cleaning techniques 
discussed in section ^ We slightly modify the cylindrical shock problem of section 15.21 
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Figure 4. This figure shows solutions of the cylindrical shock and relativistic rotor test 
problems at t = 0.4. The top two frames show P and B 2 for the cylindrical shock problem, 
and the bottom two frames show the same variables for the relativistic rotor problem. Details 
about the evolutions are given in the text. 


and monitor V • B during the subsequent evolution for (1) a free evolution (no divergence 
cleaning), (2) an evolution with hyperbolic divergence cleaning, and (3) elliptic divergence 
cleaning evolutions with D + D_ and DqDq discrete divergence operators. 

The cylindrical shock problem is modified by changing the background pressure to 
P = 1, and varying the central pressure, P c . In the examples presented here, the central 
pressure is P c = 10 and P c = 1000. We examined other values of P c , but the results were 
essentially the same as in these two cases. (The differences arise only in the overall scale of 
the constraint violation, not in the relative performance of each technique.) The evolutions are 
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Figure 5. This figure compares the hyperbolic and elliptic divergence cleaning methods in 
two different cylindrical shock problems. Each figure plots the L2 norm of V • B as a function 
of time for a free evolution, an evolution using hyperbolic divergence cleaning, and elliptic 
divergence cleaning using both D+D— and Do Do operators. Initial data for the left frame 
have a central pressure P c = 10 and the central pressure for the right frame is P c = 1000. 
Unfortunately, these L2 norms are dominated by a few points primarily near shocks, but this 
gives some indication of how violations of the constraint vary in time. 



Iflg.ol( D o) t B 1 l 



Figure 6. This figure shows the distribution of | V • B| at an instant of time, t = 0.36, for 
the cylindrical shock problem with P c = 1000. On a given grid, the number of points, N, 
with |V • B| in a specified range are counted, and log 10 (l + N) is plotted on the vertical 
axis. The total number of points is 501 2 . In the left frame the solid line corresponds to the free 
evolution, the dotted line an evolution with hyperbolic divergence cleaning, and the dashed line 
corresponds to elliptic divergence cleaning with the three-point stencil. In the right frame, the 
lines are ordered from top to bottom: free evolution, hyperbolic divergence cleaning, three- 
point elliptic divergence cleaning. Points with |V ■ B| < 10 —11 were excluded from the 
distribution on the left. 
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performed on a single grid with 501 2 points with coordinates (x,y) £ [— i, i]. During the 
evolutions V B is calculated using the central discrete operator D 0 . To emphasize the discrete 
nature of these calculations, we write the numerically calculated divergences as ( D 0 )iB l . 

Comparing the divergence cleaning techniques is more difficult than it appears at first 
blush. The first impulse is to simply plot L2 norms of (Do)iB 1 as a function of time, as shown 
in figure 0 Glancing quickly at this figure, one might conclude that the elliptic divergence 
cleaning with DoDq gives ideal results, as ||(.Do)i-B I ||2 can be made as small as desired. 
However, as discussed previously, “exact” satisfaction of the constant can be a red herring, 
in that it does not imply that the continuum constraint is “exactly” satisfied. Using the same 
elliptic projection method with a different discrete divergence operator, the I) + D_ operator, 
gives a non-zero value for 11 (/A) ),/;>' 1 12 and a better indication of error in the solution. 

A second difficulty is that the largest constraint violations appear, not surprisingly, 
near shocks. Indeed, examination of the data show that a few points near the shock 
completely dominate the values of ||(Z?o)iB*|| 2 . This makes the comparison of L2 norms 
in figure 0 problematic, as they provide almost no information about constraint violations 
in smooth parts of the solutions, where comparisons between techniques may be more 
meaningful. Moreover, since convergence in the sense of Richardson extrapolation can 
not be defined for discontinuous solutions, the norms | \(Do)iB' l \ I 2 do not become smaller 
with finer resolution, rather, the opposite occurs. With finer resolution, the shock profile is 
sharpened, and derivatives of the discontinuous variable approach the continuum derivative: 
d/dx[9[x — x')\ = 6(x — x'). Ideally one could remove points near the shock from the 
comparisons, by either tracking the shocks or simply removing points where |(_Do)i-B*| is 
judged to be too large. We have no facility for the former, and the latter strikes us as too 
arbitrary. Since L2 norms of (Do)iB l give at best limited information, in figure[6]we plot the 
distribution of \(Do)iB l \ at a single instance of time for the cylindrical shock problem with 
P c = 1000. 

6. Conclusion 

The numerical scheme presented here is for solving the relativistic MHD equations on 
multiple domains with overlapping grids. While we have not presented data from black hole 
spacetimes in this paper, it is the target application that has influenced our design decisions. 
First, we choose to work with multiple domains since excision is most naturally implemented 
with smooth boundaries adapted to the event horizon’s geometry. The flexibility of the 
overlapping grid approach allows one to easily use high resolution shock-capturing schemes 
on multiple domains. 

Secondly, we choose an ENO method with a finite difference discretization to simplify 
the transfer of information from one grid to another. While a conservative scheme for 
overlapping grids could be used E>1, working with point values instead of cell averages 
simplifies the interpolation process for arbitrary grids. ENO finite difference schemes are 
also easily extended to higher dimensions and higher orders of accuracy. 

Thirdly, we choose a central scheme to solve the RMHD equations. Central schemes are 
very efficient HRSC methods, and ideal for combining with AMR to resolve small features. 

Fourth, we investigated the use of hyperbolic and elliptic divergence cleaning to maintain 
the V • B = 0 constraint for MHD. These techniques are readily used on domains with 
arbitrary overlapping grids. We found that hyperbolic divergence cleaning often gives 
acceptable results, especially for moderate shocks. Elliptic divergence cleaning is more 
robust, but also more more computationally expensive. 


Relativistic MHD and black hole excision 


19 


Finally, it is natural to use a finite difference formulation of the fluid equations when also 
solving the finite difference Einstein equations with AMR. When both the fluid and geometric 
variables are refined in the same manner, the inconvenience of using staggered grids with 
AMR is eliminated. 
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Appendix A. CENO Reconstruction 


This appendix summarizes the CENO reconstruction scheme used by Del Zanna, Bucciantini 
and Londrillo 13011331 . based on the original scheme of Liu and Osher l [2911 . The numerical 
fluxes are constructed dimension by dimension, thus the basic algorithm is one dimensional. 
Consider an uniform grid Xi = iAx with the function i\ = v(x,). The standard one-sided 
and centered discrete differential operators are 


( D±v)i = ± 


Vi±l - Vi 


( D 0 v)i = 


Vi+l — Vi -1 


(A.l) 


Ax ’ 2A:r 

To reconstruct V{ on the interval / 2 ,2',:+]/a] we first create a linear TVD interpolating 
polynomial 

(A.2) 


= Vi + v[{x - xi), 


where v'i is the limited slope at v[ is 

v[ = minmod(U_Uj, D + Vi). (A.3) 

where the minmod limiter is defined in Other TVD limiters can be considered, such as 
the monotonized central difference limiter, but we do not consider them here. The first order 
reconstruction is v{x) = t/ 1 ' (x), which is equivalent to the TVD reconstruction, and results 
in a second-order scheme. 

Higher order reconstructions proceed hierarchically using the ENO philosophy of 
constructing multiple candidate polynomials, and then choosing the polynomial that is closest 
to the lower order polynomial. For example three candidate quadratic polynomials, Q k {x), 
A; = —1,0,1, for a second order reconstruction are 

Qi (*^) — Vi+k “h DoVi-\-k{x Xi- (_fc) T — D^-D— Vi-\-k{x Xi-\-k')“. (A.4) 

These second order polynomials are compared to the first order polynomial at the point of 
interest, x, to calculate the weighted differences 

d k (x) = a k ^Qi (x) — v^ (x)^ . (A.5) 

The weights a k are chosen to be or 1 = a 1 = 1 and a 0 = 0.7, biasing the interpolation 
towards the centered polynomial. When all d k have the same sign, the second order 
reconstruction is v(x) = Q" (x), where a is the index corresponding to the weighted 
difference with the smallest magnitude, d a (x) = min(|d; s (a:)|). When dk{x) have differing 
signs, we revert to a first order reconstruction, v(x) = v^^x). This comparison to the first 
order reconstruction yields results similar to TVD schemes near discontinuities, but gives an 
higher order reconstruction for smooth functions. 
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Figure Bl. A five point interpolation stencil is divided into three substencils, So, Si and 
S 2 , each with three points. Interpolants are calculated using all three stencils, which are then 
combined in a weighted convex sum. The nonlinear weights depend on the smoothness of the 
function. 


Appendix B. WENO interpolation 


Boundary data at interfaces between overlapping grids are obtained by interpolation. 
Following Sebastian and Shu EEZ). we have investigated both Lagrangian and WENO 
interpolation at these grid interfaces. Lagrangian interpolating polynomials work well for 
smooth functions. Near discontinuities, however, they become oscillatory, which can lead to 
unphysical states in relativistic MHD, i.e., physical primitive variables can not be obtained 
from the interpolated conservative variables. WENO interpolation avoids these oscillations 
near discontinuities by adjusting the interpolation stencil according to the local smoothness of 
the data. The WENO interpolant is constructed as a convex sum of lower order interpolations 
computed on substencils. The contribution from each substencil is weighted nonlinearly by 
measures of the function’s local smoothness. For smooth functions, WENO interpolation 
closely approximates Lagrangian interpolation. A lower-order interpolant is calculated near 
discontinuities, which is biased towards substencils where the function is smooth. In this 
appendix we summarize the WENO interpolation algorithm E3 that we use exclusively in 
our code. 

Consider a discrete function it, defined at (2 k — 1) points, Xi, i = 0, ... ,2k — 2. The 
Lagrangian interpolation polynomial ul (x) is 

2k—2 2k—2 

u L {x) = ^2 Ci{x)ui, Ci(x) = _ J . (B.l) 

2 — 0 j =0 J 

To motivate WENO interpolation, we divide the (2k — 1) points into k substencils, see 
figure IbTI and write ul(x ) as a sum of lower order interpolating polynomials on each 
substencil. Let the substencils be S r (i) = { Xi- r ,..., Xi- r +k-i} for r = 0,..., k — 1. 
The Lagrangian interpolating polynomial on each substencil is 


,( r h 


k-1 


k -1 

(x') — 'y ^ U'i— r +jC r j (x) , C rj (x) — 

j =o t=0 

We can expand ul(x ) in terms of u2(x) as 

k -1 

Ul(x) = Y dr(x)u { -[ ) (x) 


X Xi- r .\-£ 

. Xi-r-\-j Xi—r-\-t 


(B.2) 


(B.3) 


r—0 

where d r (x) are constants, or linear weights, that depend on x. Consistency requires that 

Er=0 d r ( x ) = !■ 

The Lagrangian interpolation polynomial ul(x) is written above as a sum of lower- 
order polynomials in IB. 3k WENO interpolation generalizes this by creating a convex sum. 
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Figure B2. Interpolations in two dimensions are calculated as a series of one dimensional 
interpolations. The dashed diagonal lines represent coordinate lines for a grid that requires 
interpolated at their intersection, indicated by the solid circle. The interpolation data are 
calculated from the grid with rectangular coordinate lines by first interpolating horizontally, 
using the points indicated by open circles, to obtain data at the points indicated by solid 
squares. The stencil indicated by solid squares is then used to obtain interpolation data at 
the required point. 


where the weights are nonlinearly dependent on the local smoothness of Ui. The WENO 
interpolation polynomial has the form 


k-1 

U w (x ) = U! r (x)uf / \x) 

r =0 


(B.4) 


where uj r (x) are the nonlinear weights, and for consistency we require Y2r=Q UJr ( x ) = 1- 
Following the fundamental reconstruction procedure for WENO evolution schemes, the 
weights are chosen to be 


ujr (tr) 


£j r (x) 

sr^k — 1 ~ / \ ’ 

T,s=0 Us{x) 


UJ r (x) 


dr{x) 

(e + P r {x)) 2 ' 


(B.5) 


where e is a small number that we set as e = 10 6 . The coefficients d r [x) are obtained from 
and the functions B r (x) are smoothness indicators given by 



Ax 2 '" 1 



dx, 


(B.6) 


where the limits of integration are over different substencils. All interpolations for the grids 
used in this paper are centered, x € [&*_i/ 2 , x i+ i/ 2 }. In this case the smoothness indicators 
become 1571 

Po = (10 u 2 - 31uiU i+ i + 25 u 2 +1 + lliiiUi +2 - 19u i+1 u i+2 + 4m^ +2 )/3 (B.7) 

Pi = (Bu 2 i _ l - \3ui-iiui + 13 u 2 + 5uj_iu i+ i - 13uiUi + i + 4m? +1 )/3 (B.8 ) 

P 2 = (4rti_ 2 “ 19Mi_ 2 Wi-i + 25 u 2 _ ± + llu i _ 2 Mi - SlUj-iWj + 10u?)/3 (B.9) 


Since the interpolation coefficients u> r (x) depend on the smoothness of the function, 
they must be calculated for each function individually, whereas Lagrangian interpolation 
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coefficients are only position dependent. This makes WENO interpolation more expensive 
to implement when several functions must be interpolated at a single point. Finally, 
Interpolations in two dimensions are calculated as a series of one dimensional WENO 
interpolations, as shown in figure lB2l 
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